Same as Dataset 2f but sample size of 100 and number of taxa to 100. All different interaction matrices
library(tidyverse)
── Attaching core tidyverse packages ────────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.1.4 ✔ readr 2.1.5
✔ forcats 1.0.0 ✔ stringr 1.5.1
✔ ggplot2 3.5.1 ✔ tibble 3.2.1
✔ lubridate 1.9.3 ✔ tidyr 1.3.1
✔ purrr 1.0.2 ── Conflicts ──────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors
library(igraph)
Attaching package: ‘igraph’
The following objects are masked from ‘package:lubridate’:
%--%, union
The following objects are masked from ‘package:dplyr’:
as_data_frame, groups, union
The following objects are masked from ‘package:purrr’:
compose, simplify
The following object is masked from ‘package:tidyr’:
crossing
The following object is masked from ‘package:tibble’:
as_data_frame
The following objects are masked from ‘package:stats’:
decompose, spectrum
The following object is masked from ‘package:base’:
union
library(NBZIMM)
Attaching package: ‘NBZIMM’
The following object is masked from ‘package:stringr’:
fixed
library(SpiecEasi)
Attaching package: ‘SpiecEasi’
The following object is masked from ‘package:igraph’:
make_graph
library(LIMON)
library(here)
here() starts at /Users/bealab/Desktop/Long_networks
library(lme4)
Loading required package: Matrix
Attaching package: ‘Matrix’
The following objects are masked from ‘package:SpiecEasi’:
tril, triu
The following objects are masked from ‘package:tidyr’:
expand, pack, unpack
library(Matrix)
library(tscount)
library(patchwork)
library(MASS)
Attaching package: ‘MASS’
The following object is masked from ‘package:patchwork’:
area
The following object is masked from ‘package:SpiecEasi’:
fitdistr
The following object is masked from ‘package:dplyr’:
select
library(matrixcalc)
Attaching package: ‘matrixcalc’
The following object is masked from ‘package:igraph’:
%s%
library(gridExtra)
Attaching package: ‘gridExtra’
The following object is masked from ‘package:dplyr’:
combine
library(devtools)
Loading required package: usethis
Registered S3 method overwritten by 'htmlwidgets':
method from
print.htmlwidget tools:rstudio
library(miaSim)
Loading required package: TreeSummarizedExperiment
Loading required package: SingleCellExperiment
Loading required package: SummarizedExperiment
Loading required package: MatrixGenerics
Loading required package: matrixStats
Attaching package: ‘matrixStats’
The following object is masked from ‘package:dplyr’:
count
Attaching package: ‘MatrixGenerics’
The following objects are masked from ‘package:matrixStats’:
colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads,
colMaxs, colMeans2, colMedians, colMins, colOrderStats, colProds,
colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2,
colTabulates, colVarDiffs, colVars, colWeightedMads,
colWeightedMeans, colWeightedMedians, colWeightedSds,
colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
rowWeightedSds, rowWeightedVars
Loading required package: GenomicRanges
Loading required package: stats4
Loading required package: BiocGenerics
Attaching package: ‘BiocGenerics’
The following object is masked from ‘package:gridExtra’:
combine
The following objects are masked from ‘package:igraph’:
normalize, path, union
The following objects are masked from ‘package:lubridate’:
intersect, setdiff, union
The following objects are masked from ‘package:dplyr’:
combine, intersect, setdiff, union
The following objects are masked from ‘package:stats’:
IQR, mad, sd, var, xtabs
The following objects are masked from ‘package:base’:
anyDuplicated, aperm, append, as.data.frame, basename, cbind,
colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, Position,
rank, rbind, Reduce, rownames, sapply, setdiff, table, tapply,
union, unique, unsplit, which.max, which.min
Loading required package: S4Vectors
Attaching package: ‘S4Vectors’
The following objects are masked from ‘package:Matrix’:
expand, unname
The following objects are masked from ‘package:lubridate’:
second, second<-
The following objects are masked from ‘package:dplyr’:
first, rename
The following object is masked from ‘package:tidyr’:
expand
The following object is masked from ‘package:utils’:
findMatches
The following objects are masked from ‘package:base’:
expand.grid, I, unname
Loading required package: IRanges
Attaching package: ‘IRanges’
The following object is masked from ‘package:lubridate’:
%within%
The following objects are masked from ‘package:dplyr’:
collapse, desc, slice
The following object is masked from ‘package:purrr’:
reduce
Loading required package: GenomeInfoDb
Loading required package: Biobase
Welcome to Bioconductor
Vignettes contain introductory material; view with
'browseVignettes()'. To cite Bioconductor, see
'citation("Biobase")', and for packages 'citation("pkgname")'.
Attaching package: ‘Biobase’
The following object is masked from ‘package:MatrixGenerics’:
rowMedians
The following objects are masked from ‘package:matrixStats’:
anyMissing, rowMedians
Loading required package: Biostrings
Loading required package: XVector
Attaching package: ‘XVector’
The following object is masked from ‘package:purrr’:
compact
Attaching package: ‘Biostrings’
The following object is masked from ‘package:base’:
strsplit
library(reshape2)
Attaching package: ‘reshape2’
The following object is masked from ‘package:tidyr’:
smiths
library(corpcor)
Attaching package: ‘corpcor’
The following object is masked from ‘package:matrixcalc’:
is.positive.definite
Simulate a different starting covariance matrix and counts per subject. ie each one has their own seed
# Initalize an empty list
####################################################
count_tables1 <- list()
# The Loop
####################################################
# Set the indices
for (i in 1:100) {
# Set the seed for each subject
set.seed(12345 + i)
# Create Starting covariance matrix
n <- 100^2
lower_values <- runif(n / 2, min = -0.7, max = -0.2)
upper_values <- runif(n / 2, min = 0.2, max = 0.7)
user_interactions <- c(lower_values, upper_values)
user_interactions <- sample(user_interactions)
user_A <- randomA(
n_species = 100,
interactions = user_interactions,
diagonal = -1.0,
connectance = 0.5,
scale_off_diagonal = 0.2,
symmetric = TRUE)
# Heatmap of the covariance matrix
color_palette <- colorRampPalette(c("red", "white", "blue"))(20)
heatmap(
user_A,
Colv = NA,
Rowv = NA,
col = color_palette,
zlim = c(-1, 1),
main = paste("Subject", i)
)
# Run the simulation
tse_glv <- simulateGLV(n_species = 100,
A = user_A,
t_start = 0,
t_store = 4,
stochastic = FALSE,
norm = FALSE,
error_variance = 0.1)
# Get the count table
sim_data <- tse_glv@assays@data@listData[["counts"]]
# Store the count table in the list
count_tables1[[i]] <- t(sim_data)
}
# Merge together
####################################################
# Combine the count tables
combined_count_table <- do.call(rbind, count_tables1)
# Add Rownames
rownames(combined_count_table) <- paste0("Sbj", rep(1:100,
each = nrow(count_tables1[[1]])), "_Time", 1:4)
Make Metadata and merge with the count data
# Set seed
set.seed(12345)
# Df 1 is Metadata
########################################################
meta_data <- expand.grid(Time = 1:4,ID = 1:100)
rownames(meta_data) <- rownames(combined_count_table)
# Df 2 is Metadata merged with Counts
########################################################
#Round off and increase
combined_count_table <- as.data.frame(combined_count_table + abs(min(combined_count_table)))
combined_count_table <- (combined_count_table)*10
meta_counts <- base::merge(meta_data, combined_count_table, by ="row.names",
all = TRUE)
meta_counts <- column_to_rownames(meta_counts, "Row.names")
Plot the data to look at the differences among subjects
Graphs to Check
# Individual Species Plots
########################################################
# Pivot to long data
count_long <- tidyr::pivot_longer(meta_counts, cols = starts_with("sp"), names_to = "Species")
# Plot the data
count_long %>%
ggplot(aes(x = Time, y = value, colour = as.factor(ID),
group = as.factor(ID), linetype = as.factor(ID))) +
geom_line() +
geom_point() +
geom_jitter() +
ylab("Count") +
labs(linetype = "ID", color = "ID") +
facet_wrap(~ Species) + # Create a panel for each species
theme(legend.position = "none") +
ggtitle("Time Series Data")
# Plot only two
########################################################
# Pivot to long data
metacounts_filt <- meta_counts[,c(1:6)]
count_long <- tidyr::pivot_longer(metacounts_filt, cols = starts_with("sp"), names_to = "Species")
# Plot the data
count_long %>%
ggplot(aes(x = Time, y = value, colour = as.factor(ID),
group = as.factor(ID), linetype = as.factor(ID))) +
geom_line() +
geom_point() +
geom_jitter() +
ylab("Count") +
labs(linetype = "ID", color = "ID") +
facet_wrap(~ Species) + # Create a panel for each species
theme_linedraw() +
theme(axis.text.x = element_text(family = "arial",color = "black", size = 14),
axis.text.y = element_text(family = "arial",color = "black", size = 14),
axis.title.x = element_text(family = "arial",color = "black", size = 14),
axis.title.y = element_text(family = "arial",color = "black", size = 14),
legend.position = "none",
strip.text = element_text(face = "bold", family = "arial", color = "white", size = 15))
# Distribution of counts
########################################################
hist(as.matrix(combined_count_table),
breaks = 100, main = "Distribution of GLV Data", xlab = "Counts")
Save the Counts
# Make more like counts
meta_counts[,3:102] <- round(meta_counts[,3:102])
write.csv(meta_counts, here("Data","GLV_SimData", "Dataset_2f", "GLV_100_100sp.csv"))
# Save just timepoint 2
meta_counts_filt <- meta_counts %>% filter(Time == 2)
write.csv(meta_counts_filt, here("Data","GLV_SimData", "Dataset_2f", "GLV_100_100sp_T2.csv"))